######### CANTON-LEVEL ANALYSES  ############

library(stringr)
library(data.table)
library(readstata13)
library(lfe)

library(did)
library(xtable)

library(sf)
library(ggplot2)

rm(list=ls())
gc()
graphics.off()

#Set working directory

#Load canton-level panel data
load("./data/data_canton.RData")

##### TABLE A.3 Summary statistics - cantons ########

dat <- pancant

dat[mindateref==9999, mindateref:=NA]

vars <- c("lpopintcant","lshareurbanpopintcant","lareacant","ncomscant","distpref","distspref","distbrig","ldistroads","ldistriver","ldistforest","laltmax","lrugged","lwheat","polsoc","rebel1700s","rebelcant18001806","lostcheflieu1800cant","lostcheflieu1800district","anymiss","lmenmil1801s","lmenmil1806s","lcontfoncarr1802pc")
vars <- gsub("^l","",vars)
vars <- gsub("ost","lost",vars)
vars <- rev(vars)
vars <- c("mindateref",vars)

tab <- data.frame(var=vars, mean=NA, median=NA,min=NA,max=NA,sd=NA)
for(i in 1:length(vars)){
  tab[i,-1] <- sapply(c("mean","median","min","max","sd"), do.call, args = list(dat[year==1808 & !is.na(get(vars[i])),get(vars[i])]))
}
tab <- as.data.table(tab)
tab <- t(tab) %>% as.data.table
names(tab) <- tab[1,] %>% as.character
tab <- tab[-1,]
tab <- tab[,lapply(.SD,as.numeric),.SDcols = mindateref:popintcant]

tab[,`Cadaster year`:=round(mindateref,0)]
tab[,`Canton population`:=round(popintcant,0)]
tab[,`Share urban`:=round(shareurbanpopintcant*100,2)]
tab[,`Canton area`:=round(areacant,0)]
tab[,`N communes`:=round(ncomscant,0)]
tab[,`Dist. dept capital`:=round(distpref,0)]
tab[,`Dist. arrondt capital`:=round(distspref,0)]
tab[,`Dist. brigade`:=round(distbrig,0)]
tab[,`Dist. roads`:=round(distroads,0)]
tab[,`Dist. river`:=round(distriver,0)]
tab[,`Dist. forest`:=round(distforest,0)]
tab[,`Altitude`:=round(altmax,0)]
tab[,`Ruggedness`:=round(rugged,0)]
tab[,`Wheat agricultural suitability`:=round(wheat,0)]
tab[,`Political societies (1789-1794)`:=round(polsoc,0)]
tab[,`Social conflict 1700s`:=round(rebel1700s,2)]
tab[,`Rebellions 1800-1806`:=round(rebelcant18001806,2)]
tab[,`Modified canton 1793-1800`:=round(lostcheflieu1800cant,2)]
tab[,`Modified district 1793-1800`:=round(lostcheflieu1800district,2)]
tab[,`Share men in military in arrondt in 1801`:=round(menmil1801s,2)]
tab[,`Share men in military in arrondt in 1806`:=round(menmil1806s,2)]
tab[,`Missing population 1793-1806`:=round(anymiss,2)]
tab[,`Arrondt land tax per capita in 1802`:=round(contfoncarr1802pc,2)]

out <- tab[,lapply(.SD, as.character),.SDcols=`Cadaster year`:`Arrondt land tax per capita in 1802`]
out <- t(out) %>% as.data.frame
names(out) <- c("Mean","Median","Min","Max","Std. Dev.")

N <- lapply(dat[year==1808,..vars], FUN=function(x){length(x[!is.na(x)])}) %>% unlist %>% as.vector

out <- cbind(N, out)

print(xtable(out),
      include.rownames = T, 
      include.colnames = T,
      table.placement = 't!',
      type = "latex",
      only.contents = T,
      booktabs = F,
      sanitize.text.function = identity,
      file = "tables/summary_stats_canton.tex")


## FIGURE 1: Progress of the cadaster #######

pdf("./figures/histmindateref.pdf", height=6)

h <- hist(pancant[mindateref>1806 & year==1800 & mindateref<9999,mindateref], breaks=45, plot=F)
plot(h, col=c(rep("darkgray",14), rep("lightgray",25)), main="", ylab="N cantons cadastered", xlab="year", ylim=c(0,150), xlim=c(1800,1850))
legend("topright", c("Centralized cadaster","Decentralized cadaster"), pch=22, pt.bg=c("darkgray","lightgray"))
abline(h=0)

dev.off()

#### FIGURE 2: Progress of the cadaster map ########

load("./data/admingeo1808.RData")
load("./data/maps_cadaster.RData")

coms <- st_as_sf(admingeo1808[!is.na(x) & pref==1 & !(code_dept %in% c(73, 74, 6)),], coords=c("x","y"), crs="EPSG:27572") %>% st_transform(., crs=st_crs(deptf))

deptf %>% filter(!(code_dept %in% c(73, 74, 6, 27))) -> deptf
ggplot() + geom_sf(data=cantsf[!is.na(cantsf$mindateref),], aes(fill=mindateref), alpha=.9) + scale_color_distiller(palette = "Blues", direction = -1) + 
  geom_sf(data=coms, shape=3) + labs(fill="Cadaster year", shape="Department capital") + 
  geom_sf(data=deptf, fill='transparent', lwd=.5, color='black') + theme_void() 
ggsave("./figures/map_treatment.jpg", width=7, height=7)


#### FIGURE 3: Selection of cantons during the centralized cadaster period ####

pancant[popintcant==0, popintcant:=NA]
pancant[areacant==0, areacant:=NA]

distvars <- grep("^dist",names(pancant),v=T)
pancant[,(distvars):=.SD/1000, .SDcols=distvars] #convert to km

vars <- c("shareurbanpopintcant","distroads","distforest","distriver","rugged","wheat")
for(i in 1:length(vars)){
  pancant[,paste("l",vars[i],sep=""):=log(1+get(vars[i]))]
}
vars <- c("popintcant","areacant","altmax","contfoncarr1802pc","sharemil1806cant","menmil1801s","menmil1806s")
for(i in 1:length(vars)){
  pancant[,paste("l",vars[i],sep=""):=log(get(vars[i]))]
}

pancant[,rebelcant18001806:=sum(ifelse(year %in% 1800:1806,rebelcant,0)), keyby=cantid]
pancant[,lostcheflieu1800district:=max(lostcheflieu1800district), keyby=cantid]

vars <- c("lpopintcant","lshareurbanpopintcant","lareacant","ncomscant","distpref","distspref","distbrig","ldistroads","ldistriver","ldistforest","laltmax","lrugged","lwheat","polsoc","rebel1700s","rebelcant18001806","lostcheflieu1800cant","lostcheflieu1800district","anymiss","lmenmil1801s","lmenmil1806s","lcontfoncarr1802pc")
labvars <- c("population","urban","area","N communes","dist. pref.","dist. subpref.","dist. brigade","dist. roads","dist. river","dist. forest","altitude","ruggedness","wheat","political societies 1789-1794","rebellions 1700s","rebellions 1800-6","modified canton 1793-1800","modified district 1793-1800","missing population 1793-1806","% men in military 1801","% men in military 1806","arrondt land tax per capita 1802")
summary(pancant[,..vars])
cbind(vars,labvars)
vars <- rev(vars)
labvars <- rev(labvars)

length(labvars)
length(vars)

coefplot <- function(f, labvars=vars, main="", add=F){
  
  tab <- summary(f)$coefficients %>% as.data.frame
  tab <- tab[rownames(tab)!="(Intercept)",]
  
  par(mar=c(5,12,2,2)+.1)
  
  coefs <- tab$Estimate
  ses <- tab$`Cluster s.e.` 
  ind <- 1:length(coefs)
  
  if(add==F){
    plot(coefs,ind, xlim=c(-.2,.2), yaxt='n', ylab="", xlab="", main=main, pch=16)
    abline(v=0, lty=2)
    abline(h=c(1:28),col='lightgray',lwd=.5)
    segments(coefs+1.96*ses,ind,coefs-1.96*ses,ind)
    abline(v=0, lty=2)
    #lab.vars <- labvars
    axis(2, at=ind, hadj=0, line=10,labels=labvars, outer=F,cex.axis=.8, las=1,lwd=0,lwd.ticks=0)
  }
  
  if(add==T){
    ind <- ind-.4
    points(coefs,ind, xlim=c(-.2,.2), yaxt='n', ylab="", xlab="", main=main, col='gray', pch=16)
    segments(coefs+1.96*ses,ind,coefs-1.96*ses,ind, col='gray')
  }
  
}

## Main analysis (end of period 1821 or 1815)
dat <- copy(pancant)
form <- as.formula(paste("cadcent ~ ",paste("scale(",vars,")",collapse="+",sep=""),"| code_dept | 0 | code_dept"))
models <- list(
  felm(form, data=dat[year==1821,]),
  felm(form, data=dat[year==1815,]),
  felm(form, data=dat[year==1821 & sousprefcant==0,]),
  felm(form, data=dat[year==1815 & sousprefcant==0,])
)

form <- as.formula(paste("cadcent ~ ",paste("scale(",vars,")",collapse="+",sep=""),"| 0 | 0 | code_dept"))
modelsnofe <- list(
  felm(form, data=dat[year==1821,]),
  felm(form, data=dat[year==1815,]),
  felm(form, data=dat[year==1821 & sousprefcant==0,]),
  felm(form, data=dat[year==1815 & sousprefcant==0,])
)


pdf("./figures/coefplot_selection_fenofe.pdf")

coefplot(modelsnofe[[1]], labvars=labvars, main="No department fixed effects")
coefplot(modelsnofe[[3]], labvars=labvars, add=T)
legend('topright', c("All","No admin. capital"), pch=c(16,16), col=c("black","gray"), bg='white')

coefplot(models[[1]], labvars=labvars, main="Department fixed effects")
coefplot(models[[3]], labvars=labvars, add=T)
legend('topright', c("All","No admin. capital"), pch=c(16,16), col=c("black","gray"), bg='white')

dev.off()

####### Descriptive plots non fiscal variables #######

load("./data/data_canton.RData")

dt <- copy(pancant[!is.na(mindateref) & mindateref!=9999,])
dt[,highcad:=ifelse(cantid %in% pancant[cad==1 & year==1821, cantid],1,0)]
dt[,lowcad:=ifelse(cantid %in% pancant[cad==0 & year==1821, cantid],1,0)]

### FIGURE A.10 Cadaster and public works ########

pdf("./figures/plotyear_anyadoptedbatciv.pdf", height=5)
dt[highcad==1 & anymiss==0 & year %in% 1806:1840,mean(ifelse(nadoptedbatciv>0,1,0),na.rm=T), keyby=year] %>% plot(type='l', lty=2, ylab="Any public works in canton")
dt[lowcad==1 & anymiss==0,mean(ifelse(nadoptedbatciv>0,1,0),na.rm=T), keyby=year] %>% lines()
legend("topright", lty=2:1, c("Centralized cad.","Decentralized cad."), bg='white')
dev.off()


### FIGURE A.15 Cadaster and communal land ######

pdf("./figures/plotyear_admincom.pdf", height=5)
dt[highcad==1 & anymiss==0,mean(bienscomt,na.rm=T), keyby=year] %>% plot(type='l', lty=2, ylab="Communal land activities in canton")
dt[lowcad==1 & anymiss==0,mean(bienscomt,na.rm=T), keyby=year] %>% lines()
legend("topleft", lty=2:1, c("Centralized cad.","Decentralized cad."), bg='white')
dev.off()


######## DID ###########

#Chaisemartin & D'Haultfoeuille event studies: done in STATA (see did_canton)

ploteventstudy <- function(files,mains=c(""),ylab="Estimate", xlab="", ylim=NULL){
  for(i in 1:length(files)){
    if(is.null(ylab)){ylab <- files[i]}
    add <- read.dta13(paste("./stata outputs/",files[i],".dta",sep="")) %>% as.data.table
    add <- add[add$time_to_treat %in% -5:12,]
    if(is.null(ylim)){
      min <- min(add$lb_CI_95)
      max <- max(add$up_CI_95)
      ylim <- c(min,max)
    }
    plot(add[,c("time_to_treat","point_estimate")], ylim=ylim, main=mains[i], xlab=xlab[i], ylab=ylab,  type='o')
    segments(add$time_to_treat,add$lb_CI_95,add$time_to_treat, add$up_CI_95)
    abline(h=0, v=0, lty=2)
  }
}
#### FIGURE 9 Panel A Public works (impact of the centralized cadaster) ####
#### FIGURE A.11 Public works: alternative sample and specifications ####

pdf("./figures/eventstudies_18061820_anyadoptedbatciv.pdf", height=2, width=9)
files <- c("anyadoptedbatciv_18061820")
par(mfrow=c(1,3))
xlabs <- c("lead/lags cent. cadaster")
ylim <- c(-.1,.1)
ploteventstudy(files, xlab=xlabs, ylab="Estimate", ylim=ylim)
plot.new()
plot.new()
ploteventstudy(paste(files,"_nosouspref", sep=""),  xlab=xlabs, ylab="Estimate", ylim=ylim, main="Exclude admin. capitals")
ploteventstudy(paste(files,"_souspref", sep=""), xlab=xlabs, ylab="Estimate",ylim=5*ylim, main="Only admin. capitals")
ploteventstudy(paste(files,"_deptrends", sep=""),  xlab=xlabs, ylab="Estimate",ylim=ylim, main="Department trends")
dev.off()

#### FIGURE A.12 Decentralized cadaster and large public works ####
pdf("./figures/eventstudies_post21_anyadoptedbatciv.pdf", height=5, width=7)
files <- c("anyadoptedbatciv_post21")
par(mfrow=c(2,2))
xlabs <- c("lead/lags decent. cadaster")
ylim <- c(-.08,.08)
ploteventstudy(files, xlab=xlabs, ylab="Estimate", ylim=ylim, main="All cantons")
ploteventstudy(paste(files,"_nosouspref", sep=""),  xlab=xlabs, ylab="Estimate", ylim=ylim, main="Exclude admin. capitals")
ploteventstudy(paste(files,"_souspref", sep=""), xlab=xlabs, ylab="Estimate",ylim=5*ylim, main="Only admin. capitals")
ploteventstudy(paste(files,"_deptrends", sep=""),  xlab=xlabs, ylab="Estimate",ylim=ylim, main="Department trends")
dev.off()

### FIGURE A.13 Centralized cadaster and large public works #####
pdf("./figures/eventstudies_anyadoptedbatcivlarge.pdf", height=5, width=7)
files <- c("anyadoptedbatcivlarge_18061820")
par(mfrow=c(2,2))
xlabs <- c("lead/lags decent. cadaster")
ylim <- c(-.08,.08)
ploteventstudy(files, xlab=xlabs, ylab="Estimate", ylim=ylim, main="All cantons")
ploteventstudy(paste(files,"_nosouspref", sep=""),  xlab=xlabs, ylab="Estimate", ylim=ylim, main="Exclude admin. capitals")
ploteventstudy(paste(files,"_souspref", sep=""), xlab=xlabs, ylab="Estimate",ylim=5*ylim, main="Only admin. capitals")
ploteventstudy(paste(files,"_deptrends", sep=""),  xlab=xlabs, ylab="Estimate",ylim=ylim, main="Department trends")
dev.off()

#### FIGURE 9 Panel B Communal land  ####
#### FIGURE A.16 Communal land: alternative sample and specifications ####

pdf("./figures/eventstudies_bienscomt_timetrends.pdf", height=2, width=9)

xlabs <- c("lead/lags cadaster","lead/lags cent. cadaster","lead/lags decent. cadaster")
par(mfrow=c(1,3))

ploteventstudy(c("bienscomt_all","bienscomt_central","bienscomt_loconly"),main=c("","",""), xlab=xlabs, ylim=c(-.05,.05))
ploteventstudy(paste(c("bienscomt_all","bienscomt_central","bienscomt_loconly"),"_nosouspref",sep=""),main=c("","",""),xlab=xlabs, ylim=c(-.05,.05))
ploteventstudy(paste(c("bienscomt_all","bienscomt_central","bienscomt_loconly"),"_deptrends",sep=""),main=c("","",""),xlab=xlabs, ylim=c(-.15,.15))

dev.off()


######### FIGURE A.16 (C) Callaway & Sant'Anna event studies ##########

rm(list=ls())

load("./data/data_canton.RData")

pancant <- pancant[year<1846,] #Avoid whole department missing population (Moselle)

pancant <- pancant[mindateref>1806,]

pancant[,cantidn:=.GRP, keyby=cantid]

pancant[,lpop:=log(popintcant)]
pancant[,lrugged:=log(1+rugged)]
pancant[,lwheat:=log(1+wheat)]
pancant[,ldistpref:=log(1+distpref)]
pancant[,ldistparis:=log(1+distparis)]
pancant[,ldistroads:=log(1+distroads)]
pancant[,ldistforest:=log(1+distforest)]
pancant[,lnbatciv:=log(nadoptedbatciv+1)]
pancant[,lnbank:=log(nbankint+1)]
pancant[,lshareurban:=log(1+shareurbanpopintcant)]

csdidcoefs <- function(outcome="lcontfoncint", treatvar,d=pancant, cg=c("nevertreated", "notyettreated"), adj=F){
  
  d <- d[!is.na(get(outcome)) & !is.na(get(treatvar)),]
  d[,treat:=min(ifelse(get(treatvar)==1,year,NA),na.rm=T), keyby="cantid"]
  d[is.infinite(treat) , treat:=0]
  
  if(adj==T){form <- as.formula("~ lpop+lshareurban+lrugged+lwheat+ldistparis+distpref+distspref+ldistroads+ldistforest+lnbank+rebelcant+lareacant")}
  if(adj==F){form <- as.formula("~1")}
  
  out <- att_gt(yname=outcome, 
                tname="year",
                idname="cantidn",
                gname="treat",
                control_group = cg,
                clustervars = "code_dept",
                xformla = form,
                base_period = "universal",
                data=d)
  outcoefs <- aggte(out, type="dynamic", min_e = -20, max_e=30, na.rm=T) 
  return(outcoefs)
}

plotcsdidcoefs <- function(out, col=1, main="", add=F, ylim=c(-.5,.5), xlim=c(-5,10), xlab="years pre/post"){
  out <- as.data.frame(out[c("egt","att.egt","se.egt")])
  #ylim <- ylim
  out$egt <- out$egt+1
  out <- out[out$egt %in% -5:10,]
  if(add==F){
    plot(out$egt, out$att.egt, type='o', lwd=2, col=col, ylim=ylim,  ylab="Estimate", xlab=xlab, main=main, xlim=xlim)
  }
  if(add==T){
    points(out$egt, out$att.egt, col=col, lwd=2, type='o')
  }
  segments(out$egt, out$att.egt-1.96*out$se.egt, out$egt, out$att.egt+1.96*out$se.egt, col=col)
}

cs_cantondid <- function(outcome="bienscomt"){
  out <- csdidcoefs(outcome, treatvar="cad", d=pancant[year<1848 & mindateref!=9999,], cg="notyettreated")
  out_adj <- csdidcoefs(outcome, treatvar="cad", d=pancant[year<1848 & mindateref!=9999,], cg="notyettreated", adj=T)
  outcent <- csdidcoefs(outcome, treatvar="cadcent", d=pancant[year<1848 & mindateref!=9999,])
  outcent_adj <- csdidcoefs(outcome, treatvar="cadcent", d=pancant[year<1848 & mindateref!=9999,], adj=T)
  outloc <- csdidcoefs(outcome, treatvar="cad", d=pancant[year %in% 1822:1847 & cadcent==0 & mindateref!=9999,], cg="notyettreated")
  outloc_adj <- csdidcoefs(outcome, treatvar="cad", d=pancant[year %in% 1822:1847 & cadcent==0 & mindateref!=9999,], cg="notyettreated", adj=T)
  return(list(out, out_adj, outcent, outcent_adj, outloc, outloc_adj))
}

plotres <- function(out, ylim=c(-.05,.05), xlab="", main=""){
  plotcsdidcoefs(out[[1]], ylim=ylim, xlab=xlab, main = main)
  plotcsdidcoefs(out[[2]], add=T, col=2)
  abline(h=0,  v=-1, lty=2)
}

csres_admincom <- cs_cantondid("bienscomt") 

pdf("./figures/cs_admincom.pdf", height=2, width=9)
par(mfrow=c(1,3))
plotres(csres_admincom[1:2], ylim=c(-.02,.05), xlab="leads/lags cadaster")
plotres(csres_admincom[3:4], ylim=c(-.02,.05), xlab="leads/lags cent. cadaster")
plotres(csres_admincom[5:6], ylim=c(-.02,.05), xlab="leads/lags decent. cadaster")
dev.off()


######### FIGURE A.11 bottom row Callaway & Sant'Anna event studies ##########

pancant[year %in% 1806:1820, anybatciv:=ifelse(nadoptedbatciv>0,1,0)]
out <- csdidcoefs("anybatciv", treatvar="cad", d=pancant[year<1821 & mindateref!=9999,], cg="notyettreated")
out_adj <- csdidcoefs("anybatciv", treatvar="cad", d=pancant[year<1821 & mindateref!=9999,], cg="notyettreated", adj=T)

pdf("./figures/cs_batciv.pdf", height=2, width=9)
par(mfrow=c(1,3))
plotres(list(out, out_adj), ylim=c(-.1,.1), xlab="leads/lags cadaster", main="DID: Callaway & Sant'Anna (2021)")
plot.new()
plot.new()
dev.off()


##### FIGURE A14: Communal land administration ########

load("./raw data/admingeo1808.RData")
load("./data/maps_cadaster.RData")

dt <- st_as_sf(admingeo1808[!is.na(x) & yearbiencom<1849,], coords=c("x","y"), crs="EPSG:27572") %>% st_transform(., crs=st_crs(deptf))

ggplot() + geom_sf(data=dt, aes(color=yearbiencom), size=.5) + scale_color_distiller(palette = "Reds", direction = -1) + labs(color='First year with\nshifts in communal land') + 
  geom_sf(data=deptf, alpha=0) + theme_void()
ggsave("./figures/map_admincom.jpg", width=7, height=7)


##### TABLE A5 Sample communal land #####

load("./data/data_canton.RData")

pancant[,insample:=ifelse(code_dept %in% pancant[!is.na(bienscomt),code_dept],1,0)]

pancant[mindateref>1845, mindateref:=NA]
covs <- c("log(1+wheat)","log(1+rugged)","log(1+distforest)","log(1+distroads)","log(1+distparis)","log(1+distpref)","log(1+distspref)","log(1+distbrig)","log(popintcant)","ncomscant","mindateref","newtax2")
form <- as.formula(paste("insample ~ ",paste(covs, collapse="+"),"| 0 | 0 | code_dept"))
f <- felm(form, data=pancant[year==1807,])

meandepvar <- pancant[year==1807, insample] %>% mean(.,na.rm=T)

stargazer(f, 
          covariate.labels = c("wheat suitability", "ruggedness","dist. forest","dist. roads","dist. Paris","dist. prefecture","dist. sous-prefecture","dist. gendarmes","canton population","N communes in canton","cadaster date","dept. fiscal burden"), 
          dep.var.labels = "Data on communal land activities",
          no.space = T, single.row = T, 
          omit.table.layout = "n",
          add.lines = list(c("Mean dep. var", formatC(meandepvar))),
          type="latex", float=F, out="./tables/balance_admincom.tex") 



